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To avoid instabilities in the continuum semi-classical limit of loop quantum cosmol- 
ogy models, refinement of the underlying lattice is necessary. The lattice refinement 
leads to new dynamical difference equations which, in general, do not have a uni- 
form step-size, implying complications in their analysis and solutions. We propose a 
numerical method based on Taylor expansions, which can give us the necessary in- 
formation to calculate the wave-function at any given lattice point. The method we 
propose can be applied in any lattice-refined model, while in addition the accuracy 
of the method can be estimated. Moreover, we confirm numerically the stability 
criterion which was earlier found following a von Neumann analysis. Finally, the 
'motion' of the wave-function due to the underlying discreteness of the space-time is 
investigated, for both a constant lattice, as well as lattice refinement models. 

PACS numbers: 04.60.Kz, 04.60.Pp, 98.80.Qc 



I. INTRODUCTION 

Loop quantum gravity [l| canonically quantises space-time via triads and holonomies of 
the Ashtekar connection. Whilst a full understanding of the theory has yet to be reached, 
symmetry reduced versions akin to the Wheeler-de Witt mini-superspace model have been 
successfully developed [2J. As a first approximation the quantised holonomies were taken 
to be shift operators with a fixed magnitude. This results in the quantised Hamiltonian 
constraint being a difference equation with a constant interval between points on the lattice. 
Whilst these models are reasonably successful in studying certain aspects of the quantum 
regime @, S], it has been shown that they lead to serious instabilities in the continuum, 
semi-classical limit 0, @] . In the underlying loop quantum gravity theory, the contributions 
to the (discrete) Hamiltonian operator depend on the state which describes the universe. 
As the volume grows (the universe expands), the number of contributions increases. Thus, 
the Hamiltonian constraint operator is expected to create new vertices of a lattice state 
(in addition to changing their edge labels), which in loop quantum cosmology result in a 
refinement of the discrete lattice. 

It has been recently shown how this lattice refinement effect can be modelled and how 
this approach eliminates the problematic instabilities in the continuum era Whilst the 
continuum limit of these lattice refining models can be taken, there is a complication in 
directly evolving two-dimensional wave-functions, such as those necessary to study Bianchi 
models or black hole interiors. The information needed to calculate the wave-function at 
a given lattice point is not provided by previous iterations. Recently, it has been demon- 
strated 0] that a simple local interpolation scheme can be used to approximate the necessary 
data points, allowing direct numerical evolution of such two-dimensional systems. 

In this paper, we show how Taylor expansions can be used to perform this interpolation 
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with a well-defined and predicable accuracy. We first develop the scheme for the one- 
dimensional homogeneous, isotropic cosmological case, which has analytic solutions. We 
use this simple example to show that the accuracy of the system is well controlled. We 
then study the two-dimensional case of a Schwarszchild interior, which cannot be exactly 
solved (for a general lattice refinement scheme). As our interpolation system is based only on 
Taylor expansions, all we require is that the function used to model the lattice refinement be 
analytic. This allows us to look in detail at the effects of general lattice refinement models, 
beyond the simplest cases considered thus far. We also examine the instability found in 
Ref. 0] and show numerically that the analytic conditions for stability found for a specific 
lattice refinement model, are indeed valid. We then look at how these conditions change with 
different lattice refinement models. Finally, we comment of the fact that the discrete nature 
of the underlying lattice introduces a twist into otherwise straight Gaussian wave-packets 
and explain its origin. 



II. ELEMENTS OF LOOP QUANTUM COSMOLOGY 

Loop quantum cosmology is not formulated in terms of metrics and coordinates, rather 
SU(2) holonomies of the connection, hk, and triads, p, are used. In this set up the gravita- 
tional part of the Hamiltonian constraint 1 , assuming flat homogeneous and isotropic models, 
reads |8( 

4av|*> = ^|W tr E eijk (> hjijk'h- 1 ) l*> , (2-1) 

K 7 ^ ijk 

where n = 8itG and V = \p\ 3 ^ 2 is the volume operator. We use the irreducible represen- 
tation, J = 1/2, so that the Hamiltonian constraint is not plagued of ill-behaving spurious 
solutions [9[. The action of the volume operator V on the basis states \fi) is given by 

W = bri^)=(^) 3/2 |/i). (2-2) 

The basis states are eigenstates of the triad, and eigenstates of the volume operator, 
with eigenvalues /i. They satisfy the orthonormality relation 

(/ii|/i 2 ) = <J MllW . (2.3) 

The holonomy hi (more precisely h\ ) is along the edge parallel to the basis vector, 
whose length is set by jl. The parameter 7 is the Barbero-Immirzi parameter, a constant 
ambiguity parameter that can be fixed by considering black-hole entropy calculations 
Classically, p = a 2 , where a is the usual cosmological scale factor. 



1 The reader should note that the factor ordering we are using is not the conventional one. We make 
this choice because the Hamiltonian constraint for a lattice refinement model of the form ft = /io/i~ A 
with general A, turns out to be of the same form as that for the fixed lattice case. This simplifies the 
demonstration of our numerical method to solve the quantum evolution equation of generic lattice-refined 
models. Certainly, the validity of our method is independent of the choice for the factor ordering. 
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A general state in the kinematical Hilbert space can be expanded as 

!*> = £*»• (2-4) 
A* 

with the requirement 

£V^<oo, (2.5) 
a* 

so that the state has a finite kinematical norm. 

The action of the holonomies of the Ashtekar connection on the basis states reads 



hi\n) = (csl - «0isn)|ju) , (2.6) 
where 1 is the identity 2x2 matrix, <7j are the Pauli spin matrices and cs, sh are given by 

cs|yu) = - (e^ 2 + e" 1 ^/ 2 ) = ~ (|// + /}) + |/i - //)) 

sn|/i) = l - - f^/ 2 ) \ii) = ±-.(\Li + jl)-\v-jl}) , (2.7) 

with jx a real number. 

The generalised isotropic connection |c) is defined such that, 

(c|/i) = e^ c " /2 . (2.8) 

Using Eq. (12.71) . the gravitational part of the Hamiltonian constraint, Eq. (12.1 1) , can be 
written in the form of a difference equation 

4av|*> = „ 2 l 3 ~2 E S M ~ + ^-4*] |m) , (2-9) 



where 



2K 2 /ry 3 u 2 



1 / fiyyfo 



3/x V 6 



|/i + /t| 3/2 - l/u, — /2| 3/2 . (2.10) 



From now on we set, for simplicity, 1/(2a/6k 3 ^7 5 ) = 1, so that the gravitational part of the 
Hamiltonian, Eq. (12.91) is, 

4av|^) = J2 ^ I ^ + ~^' 2 " I/* " - 2*, + * M -4/i] , (2.11) 



A 4 



and the full Hamiltonian constraint reads 

(C grwr + 7^)|*) = 0, (2.12) 
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where 7^ is the matter Hamiltonian, assumed to operate diagonally on the basis states, i.e., 
H^lfi) = Ti^ii). Using Eq. (12. lip , the full Hamiltonian constraint, Eq. fl 2 . 1 2 j) . reads 



A 3 



fi + jl\ 3/2 - |/i - /2| 3/2 [^+4^ - 2^ + ^_ 4/i ] = -H* , (2.13) 



Until recently, the magnitude of the shift operator, /i, was taken to be a constant. This 
gave the resulting difference equation a fixed step-size making its analysis much simpler. 
However, in the full theory one expects that fx will be a decreasing function of fi and hence 
the step-size in the difference equation, Eq. (I2.13p . will vary. In essence, the difference 
equation is defined on a refining lattice. It has been shown that modelling of this lattice 
refinement is crucial for the stability of the classical cosmological wave-function 

In the case of an one-dimensional system, such as the one under consideration here, the 
problem can be mapped onto a fixed lattice simply by a change of basis. This depends 
somewhat on the precise form of lattice refinement, however all that is required for it to be 
possible is that the integral f jl (/i) d[i exists. To be explicit, consider a lattice refinement 
model of the form 

—A 

where fiQ is some constant 0]. If we then make the change of variables 

/i — > v 



Mo(l-A) ' 

where A; is a constant, equal to the magnitude of the shift operator associated with these 
new coordinates, Eq. (12 . 13[) becomes 

^S(v) [V v+4k - m v + *„_ 4fe ] = -H* , (2.14) 

where 

|*> = J^» , (2.15) 

V 

and 

S(u) = \\{u + k)af 2 ' {l - A) - \(u- k)a\ 3/2/{1 - A) \ , (2.16) 

with a defined as 

/i (l- A) 
a = ^—- 

The fact that the difference equations, Eq. (12.131) and Eq. (12.141) . are of the same form is 
due to the choice of factor ordering. This allows us to directly compare numerical solutions 
of the two systems, since they will have the same large-scale limit. This is indeed the reason 
for our choice of factor ordering in Eq. (12. ip . 

Note that usually one wishes to use a self-adjoint Hamiltonian, however for simplicity 
(and to allow for a more direct comparison with the case of a black hole interior), we use the 
form given in Eq. (12. 9p . which is not self-adjoint. Certainly, the numerical method described 
below is also valid for the self-adjoint well as for different factor orderings. 
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FIG. 1: This figure demonstrates the problem of evolving a difference equation that has been 
defined on a varying lattice. Given the value of defined on two adjacent lattice points (solid 
circles) one can use the difference equation to evaluate ^ on the subsequent lattice site (open 
circle). However when one attempts to evaluate \& at the next lattice point, one finds that one 
does not have the necessary data point (square). 



III. NUMERICAL EVOLUTION OF THE DIFFERENCE EQUATION 

A. One-dimensional case 



In the previous section we have seen how to transform the Hamiltonian constraint from 
an one-dimensional difference equation defined on a varying lattice, Eq. (12.131) . to one on a 
constant lattice, Eq. (12.141) . Thus, given the initial values of \1/ on two adjacent lattice points, 
one can iterate Eq. (12. 14|) to calculate \I/ on all lattice points. However, this mapping of 
the problem onto a constant lattice is not, in general, possible in two dimensions; one needs 
to develop a system that allows the use of the difference equation defined on the varying 
lattice, Eq. (I2.13p . It is clear that, given two values of \l/ defined on two adjacent lattice 
points, one can no longer iterate the difference equation to arbitrary [i 0] (see Fig. [1] which 
illustrates this problem). Assuming that the wave-function is pre-classical [l3j], i.e., that it 
varies slowly on scales smaller than the discreteness scale, one can use a Taylor expansion 
about previously calculated lattice sites to approximate the data necessary for the next 
iteration. This is equivalent to the local interpolation method used in Ref. [7j. The main 
advantage however of the Taylor expansion method we propose here is that it allows one to 
estimate the order of the approximation and, if necessary, increase the accuracy. 

Being more explicit, consider the scheme depicted in Fig. [TJ and set 



* i , *a 



^2 , ^V+4AM 



Given the value of and ^2, one can use Eq. (I2.13P to evaluate \1>3. We then move to the 
next lattice point, so that ^2 — ^3, where the over-line indicates this is the 'new' value. To 
calculate we make a Taylor expansion about ^2 to get 



= ^ 2 + - [*3 - 2^2 + *l] [1 " A(fl)] 



+ 



\ji 0)] 2 [1 - A(»)) 



d 2 ^ 



d/i 2 



+ 



,d 3 ^ 



d/i 3 



(3.1) 



where 



AM 



AO) 



with [Mi — fj, + 4/i(/i) . 



(3.2) 



Here, we have only used the first-order terms in the expansion to calculate the required 
data point, with the second-order terms being used to keep track of the accuracy of the 
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FIG. 2: The evolution of a wave-function using the Taylor expansion scheme. The circles are the 
approximated values calculated on the varying lattice, whilst the line is the exact wave-function. 
We use (jlq = 1.0, A = —0.5, Ti^ = 1.0 x 10~ 4 , initial = 10, and initial conditions \&i = 0.0, 
# 2 = 0.5. 
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FIG. 3: The difference between the exact approximate wave- functions is less than one percent. 
This error is almost entirely due to the slight difference in the initial conditions between the two 
grids, i.e., the difference between the initial separation of and vt^- We use hq = 1.0, A = —0.5, 
TL,/, = 1.0 x 10~ 4 , initial = 1-0, and initial conditions = 0.0, ^2 = 0.5. 

approximation. This system readily extends to higher order. It should be noted that if 
higher order terms are required, then Eq. (12.131) must be used to evaluate terms like ^+8^), 
^+i2/z(/Lt) > etc., so that the higher derivatives can be calculated. 

As an example, we evaluate the wave-function using both, the exact case, Eq. (12.141) . 
and Taylor expansion method described above. A small, constant matter Hamiltonian is 
used (Hfj, = 1.0 x 10~ 4 ), to give the wave-functions some fine detail, however it should be 
remembered that there is now the possibility that gravitational back reaction may become 
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FIG. 4: The value of the next order term in the Taylor expansion. We use jiq 
Tis = 1.0 x 10 -4 , //initial = 1-0, and initial conditions = 0.0, ^2 = 0.5. 



1.0, A 



-0.5, 



important. The resulting wave-functions, obtained with both approaches, are given in Fig. El 
their difference is plotted in Fig. [31 The next order term in the Taylor expansion is shown 
in Fig. HI It is clear that, even at linear order, the Taylor expansion method is extremely 
accurate. The difference between the two wave-functions is almost entirely due to the fact 
that the separation between the two initial lattice points is not exactly the same, as a result 
of the lattice refinement of the ft (fi) scheme. This alters the initial conditions slightly, 
nevertheless the discrepancy is still less than one percent. 

It should be noted that if the lattice is not refining fast enough to ensure that the wave- 
function remains pre-classical, the interpolation method would begin to fail jfj|. However, 
this would correspond to an unstable wave-function which would not have a classical large- 
scale limit. 



B. Two-dimensional case 

The cosmological quantisation procedure used in Section [TT1 can be adapted to the 
anisotropic geometry of a black hole interior. The resulting two-dimensional Hamiltonian 
constraint is again a difference equation defined on a varying lattice 

C+ (/i, T) [\E r At +2S M ,T+2<5 T - *M-25 M ,r+25 T ] 

+C (/i, r) [(fji + 26 J ^+4fi Ml r - 2 (1 + 2 7 2 <^) /i^, r +> - 2<5 M ) ^-4^] 

5 5 2 

+ C_ (fi, T) [^^2S„r-2S T ~ ^^+25„r-2S T ] = A,r , (3.3) 

with 

C± = 26^ (y/\r±28 T \ + v^R) > (3-4) 

C = y/\ T + 6 T \-y/\T-6 T \ , (3.5) 
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FIG. 5: (a) For the fixed lattice case the two-dimensional wave-function can be calculated, given 
suitable initial conditions (solid circles), (b) In the case of a refining lattice, the data needed to 
calculate the value of the wave-function at a particular lattice site (open square) is not given by 
previous iterations (solid squares). This is due to the fact that ^(/ii,Tj) — (5 M (^_|_i, Tj) ^ 0, where 
= fj>i + 4(5 M (^j, Tj). To improve the clarity of the diagram, the explicit ,u and r dependence of 



the step-sizes is suppressed, and the notation Sr, 



-i,Tj), and S f = 5 T {ji i+ i,Ti) is used. 



where we have defined 5 M and S T as the step-sizes along the [i and r directions, respectively. 
The parameter 5, with < 5 < 1, gives the fraction of a lattice edge that the underlying 
graph changing Hamiltonian uses [5(. For clarity, the \i and r dependence in 5^ (fi,r) and 
5 T (/i, r) have been suppressed. The lattice spans the (/x, r)-plane and since fi and r are 
the coordinate lengths along the polar and radial coordinates, respectively, the (/i, r)-plane 
corresponds to the (G, r)-plane of the black-hole interior. 

Note that we have again assumed that the matter Hamiltonian acts diagonally on the 
basis states of the wave-function, namely 



(3.6) 



If Sfj, and 5 T were constant, then Eq. (I3.3P could be used to iteratively calculate the value of 
\l/ at each successive lattice point, given suitable initial conditions (see, Fig. [5^). If instead 
the lattice is refining, i.e., if 5^ and 5 T are decreasing functions of n and r, respectively, 
then we have the same problem as in the one- dimensional case (see, Fig. [5)d). As in the 
one-dimensional case, one can use Taylor expansions to calculate the necessary data points. 
In general, given a function evaluated at three (non-co-linear) coordinates, the Taylor ap- 
proximation to the value at a fourth position is given by 



/ (^4, 2/4) 



<*42 — 



dy 



%2,V2 



2&I 
dx 2 



o (si 



12) 



dy 2 



X2,V2 



(3.7) 



9 



f(xi,yi) 



/ (>4,2/4) 



°42 



°32 



7 02,2/2) 



fix < 

°12 



S X 



42 



°32 



FIG. 6: Given the value of a function evaluated at three positions, (x±,yi), (2:2,2/2) and (2:3,2/3), a 
Taylor expansion can be used to approximate the value at a fourth point, (2:4,2/4). 



where the Taylor expansion is taken about the position (£2,2/2) and we have defined 8fj = 
Xi — and <5^- = yi — yj (see, Fig. [6]). To approximate the differentials in Eq. (13. 7p . use 
points (xi,?/i) and (x 3 ,y 3 ): 



f (xi,yi) = f (x 2 ,y 2 ) + 



df 



f{x3,y 3 ) = f{x 2 ,y 2 ) +5 32 



dx 
dx 



+ 5' 



£2,2/2 



v df 
12 dy 



+ ■ 



%2,V2 



+ 51 



82 



dy 



+ ■ 



X2,V2 



(3.8) 
(3.9) 



where the dots indicate higher order terms. Solving for d x f and d y f, gives 



dx 
d£ 
dy 



X2,V2 



X2,V2 



6 -f [f (*s, 2/3) - / (x 2 , 2/2)] - ^ [/ yi ) - / (x 2 , y 2 )] + 



f If (xi,yi) - / (x 2 , y 2 )] - u -f [f (x 3 , y 3 ) - / (x 2 , y 2 )\ + 



where A = 8% 2 5 v 12 - 5 y 32 5f 2 . 

As in the one-dimensional case, higher-order terms in the Taylor expansion can be used 
to improve the accuracy of the system. Here we calculate the second-order expansion to 
demonstrate that linear interpolation is sufficient for many interesting cases. Given five 
points \l/2, ^3, ^4 and ^5, where ^ n — ^ (/x n , r n ), the second-order Taylor approximation 
for the sixth point, reads 



d$> 
dfx 



+ 



2 <9/i 2 



2 Q2, 



2 dr 2 



+ 



(3.10) 



where as before <5f- = //j — /ij and 5^ = — Tj] the derivatives are all evaluated at the point 
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FIG. 7: The initial data (solid circles) is given on two adjacent ji rows and the lattice points along 
the left-hand (small r) diagonal. This allows us to evaluate the wave-function on lattice points 
(open circles) within a triangular region of base /x max — Minitial an d height r max — Thytiai- 



(/xi,ri). We can calculate the derivatives to the necessary order by solving the following 
system of simultaneous equations: 



One gets 



dfj, 



'I' 



*1 + 



*1 + 



dfx 
dfj, 



<*31 



^5 = *1 



Ai 
A 



^ — 



A2 

A 



+ 



+ 



d 2 ^ 





2 «9 2 ^ 


+ 

l 




2 


d\i 2 


2 Or 2 i 




2 d 2 V 


+ 

l 


(S^fd 2 ^ 


2 


dy? 


2 dr 2 i 




2 d 2 V 


+ 

l 




2 


d\i 2 


2 Or 2 i 


(#i 


fd 2 ^> 


+ 


(5 5 \)V* 



<9yU 2 



<9yU 2 



A3 

A 



d 2 ^ 



Or 2 



dr 2 



A4 
A 



(3-11) 



(3.12) 



where A ; is the determinant of, 



# 1 - ^ 3 
^1 - ^ 4 

V ^1 - ^ 5 



(<^) 2 /2 (^) 2 /2\ 

(^i)V2 (5 3 T i) 2 /2 

K) 2 /2 (^Ii) 2 /2 

Ci %i (C) 2 /2 (^) 2 /2/ 



5 



11 



(3.13) 



with the (i + l) th column removed. Substituting these approximations for the differentials 
into Eq. f)3.10p . one obtains the second-order Taylor approximation to the point \1/ (/x 6 ,r 6 ), 
as required. 

This system has been implemented for the vacuum case (Ti^ = 0) for an initially Gaussian 
wave-packet along the /^-direction. As shown in Fig. [7J the initial data consist of two /z-rows 
of data, adjacent in the r-direction and the data along the left-hand diagonal. One can also 
see from Fig. [7J that only data points within a similar diagonal on the right-hand side can be 
calculated. This results in the (/i, r) being restricted to a triangle of base /i max — /^initial an d 
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FIG. 8: The wave- function is calculated by iterating the difference equation using first-order Taylor 
expansion to evaluate the necessary data points at each site. The initial conditions are a Gaussian 
on Ti n itiai = 5000, centred between //initial = 100 an d Mmax = 735, of width a = 25. The lattice 
refinement model used is t) = //~ 1//2 , 5 T (fj,,r) = r -1 / 2 . Only the relevant section of the 

wave-function is plotted and to improve clarity, only every 20 th -r lattice point is shown. In (a) the 
propagation of the the wave-function is shown, whilst in (b) the value of the second-order terms is 
plotted. The second-order corrections are typically of the order of 10 -2 % over this range. 

height r max — Tinitiai- Figure [8^ shows a typical output using the first-order approximation 
evaluated for the lattice refinement model 5^ (/i, r) = /i -1 ^ 2 , 5 T (/i, r) = r -1 / 2 , whilst Fig. [8b 
gives the second-order correction to this. It is clear that, at least for slowly varying wave- 
functions, the linear approximation is extremely accurate (higher-order corrections being 

« icr 2 %). 

IV. STABILITY OF THE SCHWARSZCHILD INTERIOR 

In Ref. |H| a von Neumann stability analysis of the difference equation, Eq. (13. 3p . was 
performed for two lattice refinement models, and it was shown that in certain circumstances 
the system is only conditionally stable. In particular, it was found that for 8^ (//, r) = /io/t -1 
and S T (/j,, t) = t t -1 , the system is unstable for jj, > It. 

We investigate this instability numerically and show that the stability condition is indeed 
correct. We do this by using the scheme described above to evaluate the wave-function, given 
two initial (consecutive) \x rows, that are empty apart from a small (10 -6 ) perturbation at 
a particular value of \x. The perturbation needs to be small to ensure we remain in the 
pre-classical regime. The system is then evaluated according to the the difference equation, 
Eq. (I3.3p . across different ranges of r. Figure fl5J) shows a typical example of how the 
amplitude of the perturbation varies with r. 

By repeating this over a range of //-positions for the perturbation, we are able to empir- 
ically confirm the stability condition \i < It. We use this method to investigate how the 
lattice refinement model alters the stability properties of the system. In particular, we find 
that for 5^ = ji~ A and 5 T = t~ a , the stability condition ranges between \i < 4r for the A = 
case (constant lattice) and \i < 1.58r for the case of A = 2.0. Figure HU1 shows the cases in 
between. In all cases, the second-order correction terms are at least an order of magnitude 
lower than the wave-functions. 

It is worth noting that for the constant lattice case, one does not have to perform any 
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FIG. 9: A perturbation of VP = 10~ 6 was put on an otherwise empty initial /x-row at /i = 230, 
for the lattice refinement model <L (a*,t) = fi~ , 5 r r) = t _1 . The maximum amplitude of the 
resulting wave-function is plotted as a function of r, for the initial perturbation set at r = 113, 
r = 114 and r = 116. Analytically, the region of stability is given for r > /i/2, i.e., the amplitude 
of the perturbation grows exponentially for r < 115 and oscillates for r > 115. Here we confirm 
this numerically. Note the different scale on the graph starting within the stable region. 
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FIG. 10: For the lattice refinement models 5^ = fx" " and 5 T = t~ , the condition of stability 
has been found. Here we plot /i/r along the critical line, i.e., the line in which the system is just 
becoming unstable, as a function of the parameter A. 
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FIG. 11: The maximum amplitude for a Gaussian wave-packet centred on /i = 1060 of width 
a = 5.25, using the lattice refinement model <5» = [Iq^/t^ 1 and b T = tqt^ 1 / 2 . For small r the 
solution does not vary significantly on scales of the order of the lattice size and hence the solution 
is stable, however as soon as the solution begins to vary this stability is lost. 

interpolations to evaluate the wave-function. Thus, one does not require the wave-function to 
be pre-classical and large changes between successive lattice points are perfectly acceptable. 
However, such wave-functions clearly do not have a semi-classical limit 2 . This is not true 
for the varying lattice cases, where for the interpolation to have a significant meaning, 
one must require that the wave-function be pre-classical (so that the derivatives in the 
Taylor expansion are small). We ensure that the stability condition holds by taking a small 
perturbation 10~ 6 ). This is demonstrated well by the fact that for the case of A = — 1, 
our results confirm the analytic considerations of Ref. [5|. 

In Ref. [s[ a second lattice refinement model was shown to be unconditionally stable, 
under certain circumstances. Specifically, for 5^ (fi, r) = /UoV^T- 4 " 1 an d S T (fi,r) = r r -1 / 2 , 
assuming the solutions do not change significantly on the scale of the step-sizes, the difference 
equation was found to be stable j5j. We find that this is indeed true, however as soon as the 
wave-functions fail to be pre-classical, i.e., as soon as there is a significant variation between 
lattice points (of the order of a few percent), they become unstable. In Fig. [HI we evaluate 
a Gaussian (centred on /i = 1060 of width a = 5.25) for this lattice refinement model. This 
solution is not entirely pre-classical, since there is a variation of ~ 0.1% between the wave- 
function evaluated on successive r lattice points. This variation grows and when it reaches 
the order of a few percent the system becomes unstable. The rate at which small initial 
variations grow, and hence the rate at which the instability becomes apparent, depends on 
the \i- and r-coordinates, growing faster for large \x and small r. 

Nevertheless, these instabilities do not represent any significant problem for loop quantum 
gravity approaches to black holes. It is already know that solutions must be pre-classical in 



2 When plotting these results, only values at the lattice points should be used, since the system is inherently 
discrete, however to guide the eye we have plotted the results with lines. 
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Region of interest 




FIG. 12: The instability present in the difference equation occurs for large n and small r. If 
one is interested in reaching the value of the wave-function at large r, one needs to have a large 
initial range of [i to have sufficient data. This can result in a portion of the domain being in the 
unstable region, which makes the code very sensitive to numerical inaccuracies. To avoid this, one 
can simply use a restricted /i-range and evaluate the wave-function as far as possible. This wave- 
function can then be used to create the initial conditions for a subsequent evaluation, provided the 
wave-function has not reached the edge of the lattice. 

order to have a well-defined continuum limit (i.e., in order for lim^oo^oo — (a*> r ) 
to be valid). What is important is that this continuum limit is always stable, which is indeed 
the case for the 5^ = noy/rn~ l , $t = r o r_1 ^ 2 lattice refinement model, but not for models of 
the form 5^ (/i, r) = fi fi' A , 5 T = (fi, r) = t t~ a . 

The presence of these instabilities leads to difficulties in the numerical implementation 
of the method described above. In particular, if one wants to calculate the wave-function 
for large values of r, then one needs to start with large range of /i, to ensure that enough 
initial data is known (see, Fig. [7]). This can mean that the system is unstable at the large 
//-side of the lattice, for small initial r. 

In practise, the value of the wave-function would usually be zero in this region, since one 
is typically interested in how an initial wave-packet evolves. Any non-zero component of 
the wave-packet in this unstable region would have met the right-hand diagonal edge of the 
lattice before the high r-region of interest is reached. However rounding errors can result 
in a non-zero perturbation, which will grow exponentially due to the inherent instability of 
the difference equation. Such difficulties can be overcome, by evaluating the wave-function 
over a small region of r. Since we now only require a smaller region of [i for this evaluation, 
the unstable region can be avoided. To reach the high r-region of interest, the system can 
be re-set using the out-putted wave-function added to a larger range of //, see Fig. (TT21) . 
Alternatively, the additional initial data necessary to reach the large value of r could be 
included at the evaluation of each /i-line, however this requires an additional three sets of 
T-data, which effectively doubles the amount of initial data required. 

V. PROPAGATION OF THE WAVE-FUNCTION 

In the case of a constant lattice, an initially centred Gaussian will move to larger //, as 
r is increased (see, Fig. (T3]). [In terms of more usual coordinates, the angle 6 increases as 
the radius increases.] The reason for this can be seen simply by considering the original 
difference equation for the simple case of a lattice point on which ^ = 1 and all other 



15 




600 800 1000 1200 600 800 1000 

(a) (b) 

FIG. 13: (a) A Gaussian centred on Qu, r) = (800, 200) twists to larger values of \i as r is increased, 
(b) As we move into the continuum limit this effect disappears. In this case the Gaussian is centred 
on (h,t) = (800,5000). Here the constant lattice case is shown however this effect persists when 
lattice refinement is modelled. 



known values are zero. Schematically, it is given by the lower hexagon in Fig. [131 centred 
around the non-zero wave-function ftT . Then, Eq. ( 13.31) with 7 = 0, implies 

^ H+25 fl ,r+2S T = 2C fi/C + . 

For r ^> 5 T , this goes to zero. To see that this implies no motion of the wave-function, 
consider the next (upper) hexagon in Fig. dH centred around the lattice point ^ / M -2<5 M ,r+25 T - 
If we 'update' the coordinates so that this central point is again called (/z, r), then, in these 
coordinates, only the lower right hand point in this (upper) hexagon is non-zero, 

In this case, Eq. (13.31) gives 

^ r A t+2<5 M ,T+2<5 T = C-/C+ , 

which tends to unity for r ^> 5 T . Thus, we find that for large r a value at one lattice point 
(yU, r), moves to one with the same /^-coordinate, (/i, r + 4<5 r ), i.e., there is no motion. 

However, when we are in a region in which the lattice discreteness is important, this is 
no longer true, as 2C /i/C + ^ and the value at one lattice point introduces a non-zero 
component to the value at a lattice point with larger fi coordinate, i.e., the wave- function 
moves to larger /i. This implies the existence of some induced rotation on the wave- function 
due to the underlying discreteness of the space-time. 

If we include lattice refinement, then the same effect occurs. Once again, there is no 
motion for r ^> 5 T and in the case of lattice refinement this requirement is reached for 
lower r, since 5 T reduces as r increases. Thus, we expect that the effect will disappear 
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FIG. 14: By considering all the lattice points on the initial row except one being zero, the 
resulting wave-function can be calculated. When the discreteness cannot be ignored, there is a 
net motion towards larger values of jx (0). The lower hexagon shows the lattice points needed to 
evaluate the first step of the wave-function, whilst the upper hexagon highlights the points needed 
for the second step. The solid circles represent lattice points at which the wave- function is unity 
whilst at the open circles the wave-function is zero. From these it can be shown that the twist 
induced on the wave-function disappears on large-scales. 



quicker than in the constant lattice case, because the lattice refinement brings us into the 
continuum limit faster. However once again, as the wave-function moves into a region in 
which the discreteness of the lattice is important, a motion will be induced. 



VI. CONCLUSIONS 

Here, we have developed a simple and intuitive prescription for evaluating two- 
dimensional wave-functions to a well-controlled level of accuracy, for arbitrary lattice re- 
finement models. We focused on black-hole interiors, however the method clearly extends 
to anisotropic Bianchi models and other systems with anisotropic symmetries. 

We have shown how the stability conditions on the Hamiltonian constraint can be inves- 
tigated using this numerical method and extended the range of lattice refinement models 
for which the stability criterion are known. 

We have also examined and explained the existence of a twist in the wave- functions, due to 
the underlying discreteness of the theory; a feature that warrants further study, particularly 
in relation to its effect in microscopic black holes. 
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